Multibubble cavitation inception 
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Center for Computational Science and E-systems, 
Japan Atomic Energy Agency, Higashi-Ueno, Taito-ku, Tokyo 110-0015, Japan 

The inception of cavitation in multibubble cases is studied numerically and theoretically to show 
that it is different from that in single-bubble cases in several aspects. Using a multibubble model 
based on the Rayleigh-Plesset equation with an acoustic interaction term, we confirmed that the 
recently reported suppression of cavitation inception due to the interaction of non-identical bubbles 
can take place not only in liquid mercury but also in water, and we found that a relatively large 
bubble can significantly decrease the cavitation threshold pressure of a nearby small bubble. By 
examining in detail the transition region where the dynamics of the suppressed bubble changes 
drastically as the inter-bubble distance changes, we determined that the explosive expansion of a 
bubble under negative pressure can be interrupted and turn into collapse even though the far-field 
liquid pressure well exceeds the bubble's threshold pressure. Numerical results suggest that the 
interruption of expansion occurs when the bubble radius is exceeded by the instantaneous unstable 
equilibrium radius of the bubble determined using the total pressure acting on the bubble. When we 
extended the discussion to systems of larger numbers of bubbles, we found that a larger number of 
bubbles have a stronger suppression effect. The present findings would be useful in understanding the 
complex behavior of cavitation bubbles in practical applications where in general many cavitation 
nuclei exist and may interact with each other. 



I. INTRODUCTION 

Cavitation in liquids is a common phenomenon that 
has been observed and used in a wide range of fields in- 
cluding mechanical and chemical engineering and nuclear 
and medical applications .— 10 1 1 1 When cavita- 
tion occurs in a liquid under negative pressure, many gas 
or vapor bubbles emerge and explosively expand. The 
cavitation bubbles then collapse violently after the neg- 
ative pressure is released, emitting high-speed liquid jets 
and/or shock waves. Because of its practical importance 
and rich physics, cavitation has long been studied after 
the pioneering work by Rayleigh. 12 However, cavitation is 
not yet fully understood, and there remain many unclear 
aspects of even its fundamental characteristics, such as 
the inception processes, threshold pressures at which cav- 
itation occurs, and lifetime of emerged bubbles i&l i 1 3 ' 1 4 1 1 5 

One of the key factors in cavitation bubble dynamics is 
bubble-bubble interaction. It has long been known that, 
when bubbles interact with each other, their dynamics 
are sometimes significantly different from that of isolated 
bubbles. In Ref. [16, for example, Smith and Mesler stud- 
ied experimentally the dynamics of a spark-generated va- 
por bubble interacting with a gas bubble in an effort to 
search for a method to reduce cavitation damage. They 
demonstrated that energy transfer between the bubbles 
alters their dynamics, and that positive pressure pulses 
emitted by the collapsing vapor bubble are reduced when 
the gas bubble stays near the vapor bubble. In Ref. 
M0rch studied the collapse of cavitation bubble clusters 
using an energy-balance equation and found that the col- 
lapse intensity depends on the volume fraction of bubbles. 
He also studied experimentally the temperature depen- 
dence of cavitation bubble dynamics and showed that 
the number of cavitation bubbles increases but their size 



decreases as the liquid temperature increases. Using a 
detailed theoretical model, Chahine and Liu 1 ^ studied 
the growth of a vapor bubble cluster in a superheated 
liquid, and found that the expansion rate of bubbles and 
the temperature drop at the bubble surface decrease as 
the number of interacting bubbles increases. In a numer- 
ical study of shock-wave propagation in a bubble cloud, 
Wang and Brennen 1 ^ showed that the expansion rate of 
bubbles at the cloud center is significantly smaller that 
that on the cloud surface. In a study of multibubble sono- 
luminescence, Mettin et al^ showed numerically that in 
a strong acoustic field, the sign of the secondary Bjerk- 
nes force (an interaction force acting between pulsating 
bubbles) depends in a complicated manner on the driv- 
ing pressure, ambient radii of bubbles, and inter-bubble 
distance. 

In this paper we examined the inception processes of 
cavitation in water in multibubble cases or, in other 
words, multi-nuclei cases. The present study was moti- 
vated by a recent numerical study of cavitation in liq- 
uid mercury and a technique to suppress it. Ida et 
a L 21 ^ 22 showed numerically that bubble-bubble interac- 
tion through pressure waves can have a significant im- 
pact on the inception of cavitation. Using a nonlinear 
multibubble model based on the Keller-Miksis equation, 
it was found that in certain situations, the explosive ex- 
pansion of a bubble (i.e., cavitation inception) is com- 
pletely suppressed by the positive pressure waves that 
larger neighboring bubbles emit when they are growing. 
This theoretical prediction was confirmed by an experi- 
ment, the results of which showed that microbubble injec- 
tion into liquid mercury suppresses cavitation inception 
and significantly reduces cavitation damage. 23 This find- 
ing implies that the inception processes of cavitation in 
multibubble cases can be significantly different from that 
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in single-bubble cases. 0.3 



We present in this paper a detailed discussion of the 
effect of bubble-bubble interaction on cavitation incep- 
tion in water. The theoretical model used in this study 
is a nonlinear multibubble model in which Rayleigh- 
Plesset equations for single bubbles are coupled through 
an interaction term that represents the bubble-emitted 
pressure waves. From a number of relevant nonlin- 
ear models^>2&2^2&2L2&22 wg gelect the c i ass i ca i 

Rayleigh-Plesset equations, which are for gas bubbles in 
an incompressible liquid, coupled through an acoustic in- 
teraction term not taking into account time delay effects 
due to the finite sound speed of water. The time history 
of the far-field liquid pressure is assumed to be a con- 
stant negative pressure following a sinusoidal decompres- 
sion. When the period of the decompression process is 
set to zero, this pressure-time history is reduced to a step 
change that has been frequently considered in cavitation 
study. Based on the model and setting, we first consider 
the interaction of two non-identical bubbles under neg- 
ative pressure. The numerical results suggest that the 
suppression of cavitation inception reported in Ref. [2l| 
can occur in water as well if the bubbles' ambient radii 
and inter-bubble distance are properly set. Then, we 
discuss several details of the results. Particular atten- 
tion is focused on the effective cavitation pressure (a dy- 
namic Blake threshold) of the suppressed bubble and the 
transition region in parameter space where the bubble's 
behavior drastically changes as, for example, the inter- 
bubblc distance changes. We also discuss the dynamics 
of larger numbers of identical bubbles to show that they 
have a stronger suppression effect. In this discussion, 
we derive an exact theoretical formula that relates the 
surface velocity of interacting bubbles to their instanta- 
neous radii. The theoretical formula is a direct extension 
of the single-bubble formula derived previously^^ and 
has high accuracy and wide applicability. Using the for- 
mula we discuss how bubble-bubble interaction changes 
the expansion rate of bubbles and the negative pressure 
amplitude in water. The obtained conclusions are con- 
sistent with previous numerical and experimental obser- 
vations. In the present investigation we found that the 
inception process of cavitation in multibubble cases can 
be much more complex than in single-bubble cases and 
that a variety of patterns of inception are possible. 



The rest of this paper is organized as follows: In Sec.HTl 
a brief review of cavitation bubble dynamics is presented, 
and in Sec. IIIII the model equations used in the present 
study are introduced. Section HVl presents numerical and 
theoretical results and discussion of these results, and 
Sec. |V] presents our concluding remarks, including com- 
ments about applications and future works. 
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FIG. 1: (Color online) Pl~R curve for Ro = 2 /im, p v = 
Pa, and k = 1. The liquid pressure pl is normalized by 
the atmospheric pressure Po. pc and Rc are the threshold 
pressure and critical radius, respectively. 



II. BRIEF REVIEW OF SINGLE-BUBBLE 
CAVITATION 

One of the well-known notions of cavitation in single- 
bubble cases (called hereafter single-bubble cavitation) 
in quasistatic cases is the Blake threshold (see, e.g., 
Refs. l2l l3lU32l ). For liquid pressures below a threshold 
value, a gas bubble (cavitation nucleus) does not have an 
equilibrium radius and thus undergoes unbounded expan- 
sion implying the occurrence of cavitation. The threshold 
liquid pressure is called the Blake threshold pressure and 
the bubble radius at the threshold is called the Blake crit- 
ical radius. These critical values are given by the pressure 
balance equation at the bubble surface, 

PL=Ph ~~R' 

where pl is the liquid pressure, pb is the internal pressure 
of the bubble given by 

P b =(Po + ^- Pv ){^-) 3 \p v , (2) 

a is the surface tension, R is the bubble radius at p^ 1 
Po = 0.1013 MPa is the atmospheric pressure, Ro the 
ambient radius of the bubble, p v the vapor pressure, and 
k the polytropic exponent of the gas inside the bubble. 
Equation ([1]) describes the balance between pl and pb 
through the surface tension force 2a/R. An example of 
the pl~R curve given by Eq. ([T|) with Rq = 2 /jm is shown 
in Fig. [TJ Here, we assumed that p v — Pa for simplicity 
and k = 1, that is, the bubble interior is isothermal. This 
result reveals that for pi, below a threshold value pc, no 
bubble radius exists that satisfies Eq. ([T]). The threshold 
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pressure, determined from dpL/dR = 0, is 



Pc = Pv 



32ct 3 



^ 27 (P + Jg-p„ )/?;; 
and the corresponding critical radius is 



exist and may interact with each other if they are suffi- 
ciently close to each other. Single-bubble study is thus 
not sufficient to fully understand cavitation dynamics in 

(3) 

practical situations. As shown in Sec. I1V1 bubble-bubble 
interaction can (sometimes drastically) alter the incep- 
tion processes of cavitation in several ways. 



Rl 



(4) 



When pl slowly decreases and exceeds pc, the bubble 
begins to explosively expand. Then, if pl holds constant 
at pl < pc, the bubble will be in a steady growth state, 
where the expansion rate, dR/dt (t being time), is nearly 
constant. The asymptotic expansion rate of such a bub- 
ble follows a simple formula*^ 



dR^ 2( Pv -p L ) 
dt ~ V 



3p 



(5) 



where p is the liquid density and, in this case, pi, is the 
liquid pressure in the far field where the bubble-emitted 
pressure is negligible. 

For pc < Pl < Pv, Eq. (TT]) has two roots which cor- 
respond to equilibrium radii.— i 30 i 31 i 32 i 33 The smaller root 
represents the stable equilibrium radius, around which 
the bubble can undergo damped oscillation. On the other 
hand, the larger root represents the unstable equilibrium 
radius, where a small deviation in bubble radius results 
in the breakdown of equilibrium: The bubble will shrink 
quickly towards the stable equilibrium radius if the de- 
viation is negative but will expand without bound if the 
deviation is positive. As shown in Sec. lIV Cl the unstable 
equilibrium radius plays an important role in our study. 

In dynamic cases where the liquid pressure varies 
rapidly and the transient motion of bubbles plays a role, 
however, the above scenario is only a rough descrip- 
tion. Detailed theoretical and numerical investigations 
of single-bubble dynamics have revealed several dynamic 
effects on cavitation inception. Researchers have found, 
for example, that even if min[pi,(£)] > pc, the instanta- 
neous bubble radius can in certain conditions be greater 
than the unstable equilibrium radius in the transient pe- 
riod, resulting in the inception of cavitation . 30 ' 31 i 33 ' 34 i 35 
This observation says that the effective threshold pres- 
sure and the effective critical radius in dynamic cases are 
different from those in the quasistatic case. Researchers 
also found that in a dynamic case the liquid viscosity, 
which obviously does not appear in the quasistatic for- 
mula, can alter the threshold pressure and the critical 
radius (see Ref. [33| and references therein). Since in the 
present study we consider the effect of the pressure waves 
emitted by growing bubbles, our problem is naturally dy- 
namic. 

As mentioned above, a number of useful insights into 
single-bubble cavitation have been published. In real- 
ity, however, not only one but many cavitation nuclei 



III. MODEL EQUATIONS 

In the present study we assume the cavitation nuclei 
to be spherical gas microbubbles. The theoretical model 
used to describe their evolution is the coupled Rayleigh- 
Plcssct equations, which read 



Ri Ri 



-Ri = -p s 
2 p 
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1 d(R)R 3 ) 
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dt 



Pb,i - -pr 



2a 



1,2,- 



4pR z 
Ri 



,N, 



~PL(t), 



(6) 



(7) 



where Ri = Ri(t) is the time-dependent radius of bubble 
i, Dij is the distance between the centers of bubbles i 
and j, pL^t) is the liquid pressure in the far field, N 
is the number of bubbles, and the overdots denote the 
time derivative d/dt. The surrounding liquid is assumed 
to be water of density p = 1000 kg/m 3 , viscosity p — 
1.002 x 10~ 3 Pa s, and surface tension a = 0.0728 N/m. 
The compressibility of water is neglected, since we are 
not interested in the details of the collapse of bubbles. 
The bubble content is assumed to be an ideal gas, and 
the internal pressure of bubble i {pb,i) is thus given by 



Pb ^ [Po + R-o [-R- 



3 k,; 



(8) 



where R^ and K>i 21X6 the ambient radius and polytropic 
exponent, respectively, of bubble i. The vapor pres- 
sure and mass exchange across the bubble surface are 
neglected. We did not consider the adiabatic behavior 
of violently collapsing bubbles, since we are interested 
only in the expansion phase. Thus, ftj is set to unity in 
most cases. Since translation of bubbles and high-order 
terms depending on the translational velocit y 36 ! 37 ' 38 ! 39 
are neglected in this model, we set the inter-bubble dis- 
tances Dij to be much greater than the ambient radii 
[e.g., as > 10(Rio + Rjo)]. Using essentially the same 
model, Bremond et al^ recently showed that the cou- 
pled Rayleigh-Plesset equation can accurately describe 
the first rapid expansion and collapse, even when the in- 
stantaneous bubble radii reach 75% of half of D^ . Based 
on this finding, we allow the maximum bubble radii to 
be comparable to their observation. 

This nonlinear system of equations describes the ra- 
dial motion of N spherical bubbles coupled through the 
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FIG. 2: (Color online) Pressure-time history assumed in the 
present study. The shown history is for T = 10 /is. A similar 
pressure profile observed experimentally can be found in, e.g., 
Ref.0. 
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FIG. 3: (Color online) Radius-time curves in single-bubble 
cases (for D12 — > 00). The ambient radii of the bubbles are 
iiio = 2 /mi and R20 = 20 /mi. The dashed line denotes the 
slope determined by Eq. ©. 



bubble-emitted pressure waves. In this system, the bub- 
bles' radial motion is driven not only by the change in 
Ph(t) but also by the pressures from the neighboring bub- 
bles described by the last term of Eq. ([5]); that is, the 
total driving pressure on bubble i (pTi) is 



JV 



PTi=PL{t)+p ^ 



1 d(R*R 



dt 



(9) 



The last term of Eq. © was derived by a well-known 
formula for the pressure wave emitted by a pulsating 
sphere, 2 ! 20 ' 21 ! 40 



P.i 



dt 



(10) 



where rj is the distance measured from the center of bub- 
ble j. 

The coupling of bubbles through pressure waves, a 
kind of bubble-bubble interaction, is known to lead 
to a variety of phenomena that single bubbles can 
never exhibit, such as attraction and repulsion of pul- 
sating bubbles^^^^i 2 .^ 4 .^ acoustic localiza- 
tion in bubbly liquids, 46,4 - avoided crossings of reso- 
nance frequencies,— large phase delays^ and filamen- 
tary structure formation in a strong sound fiekh 1 50 In 
the following section we discuss the effect of bubble- 
bubble interaction on cavitation processes. 



IV. NUMERICAL AND THEORETICAL 
INVESTIGATIONS 

A. Pressure-time history and initial conditions 

The time history of the far-field liquid pressure is as- 
sumed as follows: 



PL® = 




W{P ng 



for t < — T, 

for - T < t < 0, 

otherwise, 



with 



cos 



W 



-it- 



T) 



(11) 



(12) 



where P ng is a constant negative value and T is the pe- 
riod of the decompression process from P to P ng . This 
function represents a constant negative pressure following 
a sinusoidal decompression and is continuous up to the 
first time derivative (see Fig. [2|) . In the following discus- 
sions, P ng and T are used as control parameters. When 
one sets T — > 0, Eq. (fTTj) is reduced to a step change like 
that considered in Refs. 

0IH3O.33.35IIII 
The initial conditions assumed in the present study are 



Ri(t — — T) — R i0 , 



(13) 



Ri(t = -T)=0, (14) 
that is, the bubbles are initially at equilibrium. 

B. Competitive growth of non-identical bubbles 

Let us consider the dynamics of two non-identical bub- 
bles (bubbles 1 and 2) under negative pressure. An exam- 
ple of a single-bubble case (i.e., for D12 — > 00) is shown in 
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FIG. 4: Radius-time curves of bubble 1 for different P ng val- 
ues selected around the dynamics threshold pressure. The 
numbers shown in the panel denote P„ g /Po- 

Fig. [3] Here we set i?io = 2 /im, R 2 o — 20 /xm, Ki,2 = 1, 
P ng = — 0.25Po, and T = 10 /xs, and the corresponding 
threshold pressures are 

pci = -0.179.Po for bubble 1 (15) 

and 

p C2 = -0.007P for bubble 2. (16) 

Since P ng well exceeds the threshold pressures, both bub- 
bles undergo explosive expansion. The response of bub- 
ble 2 in the decompression process is in this case faster 
than that of bubble 1, because bubble 2 has a much 
higher threshold pressure than that of bubble 1. Hence, 
the explosive expansion of bubble 2 begins earlier than 
that of bubble 1. After the transient motion has de- 
cayed, the expansion rates of both bubbles converge to 
an almost constant value determined by Eq. ([5]). These 
observations are consistent with the well-known behav- 
ior of single cavitation bubbles. As mentioned in Sec. |TT] 
the threshold pressure in dynamic cases is in general dif- 
ferent from the value given by the quasistatic theory. 
The threshold pressure of bubble 1 in the present case 
(T = 10 (is) is slightly higher than the theoretical pre- 
diction. Figure [4] shows the radius-time curves of bubble 
1 for different P ng . The bubble cannot grow significantly 
for P ng > — 0.175Po, but it expands without bound for 
Png < — 0.176-Po- From this observation, the dynamic 
threshold pressure is found to be about — 0.176Po, which 
is only 1.7% higher than the theoretical prediction. The 
dynamic threshold comes closer to the quasistatic predic- 
tion as T increases. 

For finite -D12, the dynamics of the bubbles can have a 
different pattern. In Fig. Owe show the results for P ng = 
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FIG. 5: (Color online) Radius-time curves for P ng = 
— 0.25Po and four different D12 values. The numbers denote 
Di2/(Rio + R-2o)- The curves of R2 for different D12 values 
are indistinguishable. For D12 = 10(i?io + R20), bubble 1 can- 
not grow significantly although pl is well below the threshold 
pressure of the bubble. 



— 0.25Po with four different values of Di2- From the fig- 
ure, one finds that the expansion rate of bubble 1 is de- 
creased as D12 decreases, and the explosive expansion of 
this bubble is finally suppressed for D i2 = 10(i?io + i?2o)j 
although the negative pressure considered here clearly ex- 
ceeds its threshold pressure: the expansion ratio of bub- 
ble 1, max[i?i(t)]/i?io, for D 12 — 10(i?io + R20) is only 
about 2.57. No considerable change occurs in the dy- 
namics of bubble 2 in the shown period, because bubble 
1 is too small to cause it. This result proves that the 
suppression phenomenon reported in Refs. I2HI22I is not 
inherent in liquid mercury but can also occur in water, 
whose material properties (e.g., density, surface tension) 
are greatly different from those of mercury. In the rest 
of this subsection and in the next subsection, we discuss 
details of this phenomenon. 

The numerical result just described suggests that the 
effective threshold pressure of bubble 1 in the case of 
D12 = 10(i?io + -R20) is much lower than that predicted 
by the quasistatic theory and a more intense negative 
pressure is thus needed to cavitate bubble 1. Figure [6] 
shows the dynamics of bubble 1 for D12 = 10(i?io + R20) 
and four different values of P ng . From this, the effective 
threshold pressure of bubble 1 is deduced to be within 
the range of — 0.270Po ~ — 0.276Po, the absolute value 
of which is 1.5 times greater than the theoretical predic- 
tion. This significant change in threshold pressure is due 
to the positive pressure wave emitted by bubble 2. Bub- 
bles expanding explosively under negative pressure emit 
positive pressure waves through their radial motionji^i 
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FIG. 7: (Color online) Expansion ratios of the bubbles as 
functions of -R20. Shown are for _Rio = 2 /im, P ng — — 0.25Po, 
and t < 20 /is. 



FIG. 6: Radius-time curves of bubble 1 for D12 = 10(7?io + 
R20) and four different P ng values selected around the effec- 
tive threshold pressure of the bubble. The numbers in the 
panel denote P ng /Po- 



The positive pressure waves reduce the magnitude of the 
negative pressure in the surrounding liquid, leading to 
the need for a more intense negative (far-field) pressure 
to cavitate. The positive pressure from bubble 2 is esti- 
mated by the following simple formula: 21 



Pa(t) 



4R 2 (t) 
3£>i2 ' 



(17) 



which is given using Eqs. (| 1 0[) and (|5]) under the assump- 
tion of i?2 ~ 0. Since P ng is negative, P2(t) is positive. 
The total pressure acting on bubble 1 is thus 



P, 



P2= 1- 



4J? 2 (t) 



3D 



12 



(18) 



which is clearly higher than the far-field liquid pressure 
P ng (as shown in Sec. IIVD1 a larger number of bub- 
bles cause a greater pressure rise). Equating this to the 
threshold pressure of bubble 1 (pci), we have an ap- 
proximate formula for the effective threshold pressure of 
bubble 1: 



p — 



PCI 



4R2 



(19) 



Since the denominator of the right-hand side is smaller 
than unity, this gives a pressure lower than pci- This 
equation suggests that because the problem is essentially 
dynamic, the effective threshold pressure of a bubble can- 
not be determined uniquely but depends on the instan- 
taneous radius of the neighboring bubble. Rearranging 



Eq. (|19p , one obtains a formula for the instantaneous ra- 
dius as follows: 



(20) 



For pci = -0.179F and D 12 = IQ{R W + R20) with 
P ng = -0.25P , this formula gives R 2 = 0.213L>i 2 = 
46.9/xm. This criterion appears to be fulfilled in the ex- 
ample of Fig. [5] the radius of bubble 2 becomes greater 
than this criterion before bubble 1 begins to expand ex- 
plosively. 

The observed suppression of explosive expansion, or 
competitive growth of bubbles, occurs also for other 
couples of bubbles. In Fig. [7] we show the expansion 
ratios of bubbles, m&x[Rj(t)]/ Rjq, for i?io = 2 /jm, 
P ng — — 0.25Po7 and t < 20 /is as functions of i?20- 
Here the inter-bubble distance was fixed as D12 = 220 
/im, which corresponds to 10(i?io + -R20) in the previ- 
ous example. From this figure one finds that the ex- 
plosive expansion of bubble 1 is completely suppressed 
when R20 > 17 /im. This figure also suggests that when 
i?20 < Rio the expansion of bubble 2 can be suppressed 
by bubble 1, and that for R20 > Rw the expansion ra- 
tio of bubble 2 decreases monotonically as P20 increases. 
The latter is due to the fact that the expansion rate of 
bubble 2 has an almost constant value determined by 
Eq. |5|) although R 20 changes, and hence a larger i? 2 o 
gives a smaller expansion ratio. In Fig.[8l we show phase 
diagrams of max[i?i(i)]/i?io in a wider parameter range. 
Here, both R 2 q and D12 are used as parameters, and 
three different values of Riq [(a) 1.6, (b) 2, and (c) 3 /im] 
are assumed. The white regions seen in the right side 
of the panels denote cases where the explosive expansion 
of bubble 1 is completely suppressed. As can be clearly 
seen, the white region becomes narrower as -Rio increases. 
This is because pci increases and consequently the pres- 
sure from bubble 2 needed for the suppression of bubble 1 
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FIG. 8: (Color online) Phase diagrams of max[_Ri(t)]/7lio as 
functions of R20 and D12 for three different Rio values. 



becomes greater as Rio increases. From the same figure, 
one can also see that larger R20 and smaller D12 values 
are more effective for cavitation suppression. This is con- 
sistent with the result predicted by Eq. ([T7|) or (flU)) . and 
also with the previous numerical finding.— The above 
results prove that the suppression of explosive expansion 
can occur if the bubbles' ambient radii and inter-bubble 
distance are appropriately set. 



FIG. 9: Radius-time curves of bubble 1 for P ng = — 0.25Po 
and five different D12 values. The numbers denote Z?i2/(-Rio + 
R20). ForDi2 = 12(^10 + ^20) ~ 12.1(7210 +-R20), the bubble 
collapses although pl (the lower panel) holds constant at a 
negative value. 



C. Interrupted expansion in systems of 
non-identical bubbles 

The results shown in Fig. [5] imply that a transition of 
bubble dynamics takes place in a parameter range be- 
tween D12 = 10(i?io + -R20) and 15(i?io + #20)- Here we 
clarify what occurs in the transition region. In Fig. [9J 
we show the radius-time curves of bubble 1 for i?io = 2 
(xm, R20 = 20 fxm, P ng = — 0.25Po, an d five different 
D12 values selected from the above-mentioned param- 
eter range. The other parameters were set as in the 
above examples. As shown previously, decreasing D12 
results in the decrease of the expansion rate of bubble 
1. In the parameter range considered here, however, one 
more interacting change can be found: The expansion of 
bubble 1 is interrupted at a moment. In the results for 
D 12 = 12(R W + R 20 ) - 12.1(i?i + R20), one finds that 
bubble 1 first expands considerably, but then turns into 
collapse although holds constant at a negative value 
that exceeds the quasistatic threshold pressure. Such a 
behavior is not allowed for isolated bubbles, which can 
only expand without bound. This observation suggests 
that bubble-bubble interaction sometimes causes a sig- 
nificant change in the lifetime of cavitating bubbles. 

Let us consider the mechanism underlying this obser- 
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FIG. 10: (Color online) Total driving pressure on bubble 1 
for D12 ~ 12.05(i?io + i?2o). The dashed line denotes the 
quasistatic threshold pressure of bubble 1 (pci = — 0.179Po). 



vation. As known from Eq. (fTT|) , the total pressure acting 
on bubble 1 increases as time goes on [i.e., as R2(t) be- 
comes greater]. Indeed, the total pressure at bubble l's 
position determined numerically by 



Pti = PL 



P d(RjR 2 ) 



D 



12 



(It 



(21) 



increases for t > (see Fig. fTO]) . From this, one may pre- 
sume that the interruption of bubble expansion occurs 
when the total pressure rises above the quasistatic thresh- 
old pressure of bubble 1. This conjecture is, however, in- 
correct or insufficient. The total pressure is clearly higher 
than the threshold pressure in most periods except for a 
short duration around t = 0, and their crossing points 
are far from the time of interruption (see Fig. [T0|) . The 
notion of threshold pressure is therefore useless for the 
present purpose. 

We suggest using the unstable equilibrium radius to 
understand the interrupted expansion. As mentioned in 
Sec. [Hj a bubble has an unstable equilibrium radius un- 
der the condition of pc < Pl < Pv and begins unbounded 
expansion when the bubble's instantaneous radius be- 
comes greater than the unstable equilibrium radius, even 
though pl > pc- We reveal here that the unstable equi- 
librium radius plays an essential role in the occurrence of 
interrupted expansion. Figure [TT] shows the radius-time 
curves of bubble 1 for three different D12 values and the 
corresponding unstable equilibrium radius (called here- 
after Ruel)- Here Rjjei was determined using Eq. ([T]) 
by replacing p^ with pti (I21[) . Since pri is time depen- 
dent, Rjjei varies in time. For D12 = 12.5(i?io + i?2o)j 
Ri well exceeds Ruei for /is < t < 37 /is, and hence 
bubble 1 can expand rapidly as a single bubble does. For 
D12 = 12.05(i?io + R20), Ri is slightly larger than Ru e i 
for /is < t < 22 /is and bubble 1 expands mildly during 
this period. However, Ri is exceeded by Rjjei at about 
t = 22 /is, and then bubble 1 stops expanding and be- 
gins collapsing. For D 12 = 12(-Rio + R20), Ri is smaller 
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FIG. 11: (Color online) Radius and unstable equilibrium ra- 
dius of bubble 1 for different D12 values as functions of time. 
The thick curves denote Ri and the thin curves with circles 
denote Ruei- The numbers denote Di2/(i?io + R2o)- 



than or almost equal to Rjjei in most periods, and hence 
bubble 1 cannot expand considerably and behaves as a 
stable bubble under positive absolute pressure. As can 
be seen in Fig. [TTT b). the crossing point of R\ with Rjjei 
correlates well with the timing for bubble 1 to stop accel- 
erating. This observation proves that the instantaneous 
unstable equilibrium radius can be used as a probe for 
the interruption of bubble expansion. 

Until now we have only considered cases of T = 10 /xs. 
This time length is greater than the characteristic time 
of bubble 1, Tn = 0.296 /is for i?i = 2 /im, given by 



T, 



61 
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(22) 
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FIG. 12: (Color online) Same as Fig. El but for T = 2 /Lts and 
smaller D 12 /(Rio + 7?2o). 



with lo\ being the isothermal eigcnfrcqucncy of bubble 1: 

m - \lik( m+ £)- <23) 

Here we briefly examine the effect of smaller T. For 
smaller T, exceeds the threshold pressures of the 

bubbles earlier than for larger T and the bubbles' dynam- 
ics and interaction effect would thus be altered. Also, as 
T approaches T&i , the response of bubble 1 to the decom- 
pression process should change. Figures [T2l and fT3l show 
results for T = 2 /Lis and 0.3 /lis (~ Thi), respectively, 
with P ng = — 0.25Po- In the case of T = 2 /is, one finds 
that the inter-bubble distance of D12 = 12(i?io + R20), 
which was sufficient to cause interrupted expansion when 
T = 10 fis, is now not sufficient and a smaller _D 12 is re- 
quired. This is because the rapid expansion of bubble 




Time [|j,s] 

FIG. 13: (Color online) Same as Fig. \U\ but for T = 0.3 fis 
and smaller Di2/(Rio + R2o)- 



1 begins in an earlier stage than in the previous exam- 
ple, where the radius and expansion rate of bubble 2 are 
not so large. This result reveals that the dynamics of 
interacting bubbles depends on the time history of the 
motion. The same tendency can be found in the case of 
T = 0.3 fis, where an even smaller D12 is needed to inter- 
rupt the expansion of bubble 1. However, these results 
are qualitatively the same as that for T = 10 /is: Bubble 
1 can grow considerably only when R\ > Rjjei- 



D. The case of larger numbers of bubbles 

So far we have considered double-bubble cases. In re- 
alistic cavitation, however, a large number of bubbles 
emerge at the same time and should interact with each 
other in a complicated manner. As a first step toward 
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understanding the inception processes in many-bubble 
cases, in this subsection we discuss how a larger number 
of gas bubbles change the negative pressure amplitude in 
water. As shown below, a larger number of bubbles de- 
crease the negative pressure amplitude more significantly 
and thus have a stronger suppression effect. 

In this discussion, we categorize the bubbles consid- 
ered into two groups. The bubbles in the first group are 
assumed to be very small, and their rapid expansion is 
completely suppressed by other bubbles. The bubbles 
in the second group are relatively large and can expand 
rapidly under an assumed negative pressure. This as- 
sumption allows us to neglect the impact from the first 
group on the second group, because the amplitudes of 
the pressure waves from the bubbles in the first group 
are negligibly small. We furthermore assume, as done 
in the local homogeneity assumption,— that the bubbles 
in the second group are identical and have the same dy- 
namics. Under these assumptions, we derive a theoretical 
formula for R of the bubbles in the second group and use 
it to examine the influence of bubble-bubble interaction 
on the expansion rate and negative pressure amplitude. 
Since the expansion rate depends on the total negative 
pressure on the bubble, one can see from the theoretical 
formula how the bubbles change the negative pressure 
value. The theoretical result is verified numerically for 
few-bubble systems using a test problem similar to that 
considered by Chahine and LiuJ^ 

We derive the theoretical formula for R by integrat- 
ing the coupled Rayleigh-Plesset equation. It is known 
that the Rayleigh-Plesset equation for single bubbles can 
be integrated analytically if the external pressure is con- 
stant and the liquid viscosity is neglected ^ 30 ' 35 The de- 
rived formulas, which describe the relation between R(t) 
and R(t), have provided several important insights into 
bubble dynamics, such as the asymptotic expansion rate 
under constant negative pressure and the minimum size 
of a bubble at collapse. Here we extend the formulas 
to multibubble systems. A theoretical formula for multi- 
bubble systems was given in Ref.0 under the assumptions 
of R ss and Ri ^> Rio- The present formula is derived 
without these assumptions and thus has higher accuracy 
and wider applicability. We first consider a system of two 
identical bubbles and then extend to systems of larger 
numbers of identical bubbles arranged symmetrically. 

When N — 2 and Rx(t) = R 2 {t), Eq. © is reduced to 
a single differential equation, 

RlRl + _ El _ _ pail _ (24) 

which can be rewritten as 

^(2^1^! + 3RD + 7 ^-(2R 1 R 1 + ARl) = -p s<1 . (25) 
2 2L>i2 p 

We attempt here to transform the left-hand side of this 
equation to a time derivative. Using the identity of 
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{R" l R() = 2R X R X + mR(, (26) 



where m is an integer, the left-hand side of Eq. ([25)) is 
rewritten as 
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From Eqs. (J2TJ) and (J25J) we have 
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Assuming dph/dt — and p, = 0, the right-hand side of 
this equation can also be rewritten into a time derivative, 
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dt 

where 
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Substituting Eq. 
finally have 



log-Ri 



into Eq. 



for K% ^ 1, 
for K\ = 1. 



(29) 



(30) 



and integrating, we 



1 + -M i? = + - (Po + — F 31 

F>i2 J ip P\ Rw, 

2a a 
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Here, a is the constant of integration which is determined 
by an initial condition. For R\(t = 0) = Rio and R\{t = 
0) = 0, for instance, 
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(32) 
(33) 



for k — 1. We should note that the vapor pressure p v is 
neglected in this theory but it can easily be taken into 
account by only replacing Po-, Pl with P — p v , — p v , 
respectively. 

To confirm the accuracy of the presented formula, we 
performed a numerical test. In Fig.[TJ]we show numerical 
and theoretical results for Rio = 20 /xm, D\2 — 20i?io, 
Ki = 1, P ng = -0.25P , M = Pa s, and T — [is. The 
numerical results were obtained by directly solving the 
coupled Rayleigh-Plesset equation (j2"4"|) , and the theoret- 
ical results were obtained using Eq. (f3"Tj) with R\ given 
numerically and a at t — p,s. These results are indis- 
tinguishable, proving the accuracy of Eq. (|3"Tj) . Figure [Ml 
also shows the result for Dyi —> oo, from which one finds 
that bubble-bubble interaction decreases the expansion 
rate of the bubbles. In Fig. [T5] we show the numerical 
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FIG. 14: (Color online) Expansion rates of coupled bubbles 
for Rio — 20 fim, D\2 = 20i?io, and fi — Pa s. The lower 
panel shows pl assumed in this example, normalized by Po 
(Png = — 0.25Po, T — /Us). The numerical (the solid curve) 
and theoretical (the dots) results of the expansion rate are 
indistinguishable. The dash-dotted curve shown in the upper 
panel is the expansion rate of a single bubble with the same 
ambient radius. 




Time [|as] 



FIG. 15: (Color online) Same as Fig. HI but for fi = 1.002 x 
10 -3 Pa s. The lower panel shows the instantaneous Reynolds 
number given by Eq. (|34|l . 



result of a viscous case with // = 1.002 x 10~ 3 Pa s. 
The result is in close agreement with the theory. This 
is because the viscosity has, in the present case, only a 
minor effect in most periods; see Fig. 115( b) which reveals 
that the instantaneous Reynolds number at the bubble 
surface, defined here by 



Re 



inertial term 
viscous term 



pRi 



3 pRiRi 
8 



(34) 



is large except for in the early stage. 

In Eq. (|3"Tj) . the effect of bubble-bubble interaction is 
concentrated in the term of 



Ri 



(35) 



For D±2 — > oo, this converges to unity and the for- 
mula is reduced to that for single bubbles ^ 30 ' 35 Obvi- 
ously a smaller D12 leads to a smaller expansion rate. 
This indicates that bubble-bubble interaction decreases 
the expansion rate as observed in the above numerical 
test. This conclusion is consistent with earlier numerical 
observations /' 14 ! 15 ' 18 ^ 9 The decrease in expansion rate 
obviously indicates the reduction of the negative pres- 
sure amplitude at the bubbles' position. As discussed in 
Sec. IIVB1 this reduction is caused by the positive pres- 
sure waves emitted by the bubbles themselves^ 1 - that is, 
bubbles expanding under negative pressure emit positive 
pressure waves that reduce the negative pressure ampli- 
tude in the liquid. 

Next, we briefly consider cases of larger numbers of 
bubbles using the theoretical formula. Here, following 
Ref. [H, we discuss systems of identical bubbles arranged 
symmetrically. This configuration allows us to use a 
model equation similar to Eq. (j24|) . Examples considered 
are three bubbles arranged as a regular triangle (Case 
3), four bubbles arranged as a regular tetragon (Case 
4a), four bubbles arranged as a regular tetrahedron (Case 
4b), and eight bubbles arranged as a regular hexahedron 
(Case 8). We assume that the center-to-center distances 
between nearest-neighbor bubbles (i.e., the length of the 
sides of the regular arrangements) in all cases are the 
same, D 12 - These assumptions reduce Eq. (j6]) to a single 
differential equation, 
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dt 
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where E a is the effective inter-bubble distance given as 
follows: 
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FIG. 16: (Color online) Bubble radii and total driving pres- 
sures on bubble 1 for different bubble populations. The pa- 
rameters are the same as those used in Fig. 1141 The numbers 
shown in the panels indicate the case numbers ("1" and "2" 
are for the single- and double-bubble cases, respectively). As 
indicated by the arrows, the radius decreases and the total 
driving pressure increases as the bubble population increases. 

Since Eq. ([35]) has the same form as the double-bubble 
formula (j2~4")l , the exact formula ([51]) can be used for the 
present cases by replacing D\2 with E a . From Eq. (|37p . 
one easily finds 

D 12 >E 3 > E 4a > E 4b > E s , (38) 

that is, a larger number of bubbles or a denser popula- 
tion leads to a stronger interaction effect and a smaller 
expansion rate (note that Case 4b, whose mean inter- 
bubble distance is Z?i2, is denser than Case 4a). This 
is confirmed by the numerical result shown in Fig. 116( a) 
given by directly solving Eq. (|36p . which reveals the same 
tendency as expected. The present theoretical result is 
consistent with the numerical finding by Chahine and 
Liu^ 

The above theoretical result can also be interpreted as 
indicating that the presence of a larger number of bub- 
bles causes a greater reduction in the negative pressure 
amplitude. In a system of a larger number of bubbles, the 



total amplitude of the bubble-emitted positive pressures 
should be greater and hence the negative pressure in the 
liquid, the driving force on the bubbles, is more reduced 
compared to the case of two bubbles [see Fig. RTOT b)] . This 
naturally results in a greater reduction of expansion rate, 
as seen above. From this argument one can conclude that 
a larger number of bubbles have a stronger suppression 
effect on neighboring (smaller) bubbles than that of sin- 
gle bubbles. 

Although we have only considered few-bubble systems, 
the theoretical formula for R (f3"Tj) would be applicable to 
a larger cluster of bubbles. In recent years, there has been 
an effort to use a very simple model to study the dynam- 
ics of large bubble clusters . 14 ' 15 In those works, the degree 
of freedom of a bubble cluster is significantly reduced by 
a local homogeneity assumption— and a model equation 
very similar to Eq. (|24|) is derived. The present theo- 
retical formula provides a solution of the simple model 
equation. 



V. CONCLUDING REMARKS 

We have studied the inception processes of cavitation 
in multibubble cases, where multiple cavitation nuclei 
interact with each other, and have shown that bubble- 
bubble interaction changes the inception processes in a 
variety of ways. Performing numerical simulations of the 
dynamics of non-identical bubbles under negative pres- 
sure, we have demonstrated that the suppression of the 
explosive expansion of small bubbles by bubbles expand- 
ing earlier, recently reported for liquid mercury^ is pos- 
sible in water as well. To more deeply understand the 
numerical observation, we have discussed the effective 
threshold pressure of interacting bubbles. We found that 
even a bubble can significantly decrease (by about 50%) 
the effective threshold pressure of a smaller neighboring 
bubble. This change of threshold value is much more sig- 
nificant than that caused by the dynamic effect due to 
rapid change in the far-field liquid pressure, which is only 
about 1.7% in our case. 

From a detailed analysis of the transition region 
where the dynamics of the suppressed bubble drastically 
changes as the inter-bubble distance changes, we have 
revealed that the explosive expansion of a bubble under 
negative pressure can be interrupted and turn into col- 
lapse even though the far-field liquid pressure remains 
well below the bubble's (quasistatic) threshold pressure. 
Using the notion of unstable equilibrium radius, we have 
found that the interruption of bubble expansion takes 
place when the instantaneous bubble radius is exceeded 
by the instantaneous unstable equilibrium radius deter- 
mined using the total pressure acting on the bubble. 
Both the suppression and interruption of bubble expan- 
sion are caused by the pressure wave that a neighboring 
bubble emits when it grows. 

By analytically integrating the acoustically coupled 
Rayleigh-Plesset equations under the assumption of con- 
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stant far-field liquid pressure, we have derived an exact 
formula for the surface velocity of interacting identical 
bubbles. The formula shows clearly that bubble-bubble 
interaction decreases the expansion rate of bubbles and 
that the decrease is more prominent for systems of a 
larger number of bubbles. This theoretical result is con- 
sistent with earlier numerical findings. 18 This change in 
the expansion rate implies that the presence of a larger 
number of bubbles reduces the negative pressure ampli- 
tude in water to a greater degree than in the case with 
two bubbles, and hence a larger number of bubbles have 
a stronger suppression effect. The presented exact for- 
mula would be applicable to a large bubble cluster by 
incorporating the local homogeneity assumption. 52 

The present findings could be a key to understanding 
the complex behavior of cavitation bubbles in practical 
situations where a large number of cavitation nuclei ex- 
ist and interact with each other. From the present re- 
sults, one can for example speculate as follows: Even if 
a liquid involves a large number of cavitation nuclei and 
one imposes a strong negative pressure that can cavitate 
all of the individual nuclei, only part of them will grow 
and be observed because bubbles growing earlier sup- 
press the expansion of the remaining neighboring nuclei. 
If this scenario is true, prediction of the actual population 
of cavitation bubbles in practical situations will require 
that, among other things, one know not only the initial 
distribution of cavitation nuclei but also their dynamic 
behavior in which bubble-bubble interaction plays a sig- 
nificant role. Also, the present findings may be relevant 
to the study of the supper stability of nanobubbles, 53 where 



the interaction between nanobubbles on a hydrophobic 
surface and cavitation microbubbles is seen. In that 
study, Borkent et al. found experimentally that surface 
nanobubbles do not cavitate even for a sufficiently strong 
negative pressure, while cavitation bubbles originating 
from microscopic cracks are rapidly expanding. 

More detailed analysis of multibubble cavitation dy- 
namics by bifurcation and perturbation theories or 
other theories that have been used for single-bubble 
stufh— — — — — would be an interesting subject to pur- 
sue. Also, detailed parametric study of the effective 
threshold pressure in multibubble cases could provide 
useful insights into the cavitation threshold pressure in 
practical applications. Statistical analysis of the ac- 
tual population of explosively expanding bubbles in cases 
where non-identical nuclei interact would also be a mean- 
ingful area of study. Combining the presented pictures 
of cavitation inception with a scenario of nuclei merging, 
like that proposed by Marschall et al., 13 may provide a 
more realistic picture of multibubble cavitation incep- 
tion. 
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